
options(rasterTmpDir='C:/temp')
elec_aicd_dir <- file.path(wkdir,'local','gis_data','019_utilitiesCommunication','elec')

##
## IMPORTANT PATH to Rtools to run Cpp
##
Sys.setenv(PATH = paste("C://WBG//Utility//Rtools//bin", Sys.getenv("PATH"), sep=";"))
Sys.setenv(BINPREF = "C://WBG//Utility//Rtools//mingw_$(WIN)//bin//")

require(Rcpp)
setwd(wkdir)
output_version <- '2021-03-23'

ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01',crs=4326)

ntl.zones.prj <- ntl.zones %>%
  # 102022
  st_transform(crs="+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs") %>% 
  st_centroid() %>% 
  st_geometry()

M.ntl.zones <- ntl.zones.prj %>% 
  st_transform(4326) %>%
  st_coordinates() 

# convert polygons to points 
ntl.zones.pts = st_as_sf(as.data.frame(M.ntl.zones), coords = c("X", "Y"), crs = 4326)
ntl.zones_df = st_drop_geometry(subset(ntl.zones,select=c(OBJECTID, GID_0)))
ntl.zones.pts = dplyr::bind_cols(ntl.zones.pts,ntl.zones_df)

nCpu = 6
##
## AICD electricity
##
elec_aicd_url <- 'https://www.infrastructureafrica.org/library/doc/1045/aicd-all-countries-electricity-transmission-network/'

elec_layer <- 'AICD_ALL Countries Electricity Transmission Network'
elec_aicd <- read_sf(dsn=elec_aicd_dir,layer=elec_layer)
elec_aicd = subset(elec_aicd,AICDstatus=="Existing")
elec_aicd = elec_aicd[!st_is_empty(elec_aicd),,drop=FALSE]
ntl.zones.pts$dist2elecafr<- NA
ntl.zones.pts$dist2elecafr <- unlist(nngeo::st_nn(ntl.zones.pts, elec_aicd, k = 1, returnDist = T, parallel = nCpu)$dist) / 1000
summary(ntl.zones.pts$dist2elecafr)
elec_reg <- subset(elec_aicd,COUNTRY %in% c('CMR','NER','NGA','TCD'))
unique(elec_reg$COUNTRY)

##
## nearest in region
##
ntl.zones.pts$dist2elecreg <- unlist(nngeo::st_nn(ntl.zones.pts, elec_reg, k = 1, returnDist = T, parallel = nCpu)$dist) / 1000
summary(ntl.zones.pts$dist2elecreg)

##
## nearest in country only
##
unique(ntl.zones.pts$GID_0)
elec_cmr <- subset(elec_reg,COUNTRY=="CMR")
ntl.zones.pts$dist2elecctr[ntl.zones.pts$GID_0=="CMR"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="CMR",], 
                                                                              elec_cmr, k = 1, returnDist = T, parallel = nCpu)$dist) / 1000

elec_ner <- subset(elec_reg,COUNTRY=="NER")
ntl.zones.pts$dist2elecctr[ntl.zones.pts$GID_0=="NER"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="NER",], 
                                                                              elec_ner, k = 1, returnDist = T, parallel = nCpu)$dist) / 1000

elec_nga <- subset(elec_reg,COUNTRY=="NGA")
ntl.zones.pts$dist2elecctr[ntl.zones.pts$GID_0=="NGA"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="NGA",], 
                                                                              elec_nga, k = 1, returnDist = T, parallel = nCpu)$dist) / 1000

elec_tcd <- subset(elec_reg,COUNTRY=="TCD")
ntl.zones.pts$dist2elecctr[ntl.zones.pts$GID_0=="TCD"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="TCD",], 
                                                                              elec_tcd, k = 1, returnDist = T, parallel = nCpu)$dist) / 1000

names(ntl.zones.pts)
summary(ntl.zones.pts)

## output
save.dta13(st_drop_geometry(ntl.zones.pts), paste0(wkdir,"/proc_data/ELEC_AICD",output_version,".dta"))

